#if !defined(glib_g2_smatrix_h)
#define glib_g2_smatrix_h

#include "vector.h"

namespace g2 {

//symmetric matrix
template <class T>
class basic_smatrix;
template <class T>
class basic_matrix;

typedef basic_smatrix<double> smatrix;
typedef basic_smatrix<float> smatrixf;
typedef basic_smatrix<long double> smatrixl;

template <class T>
class basic_smatrix {
public:
    T xx, xy, yy;
public:
	typedef T * iterator;
	typedef const T * const_iterator;
public:
	T * begin()
		{ return &xx; }
	const T * begin() const
		{ return &xx; }
	T * end()
		{ return &yy + 1; }
	const T * end() const
		{ return &yy + 1; }
	static int size() 
		{ return 3; }
	static int const diagonal_elements = 2;
public:
	basic_smatrix()
		{ }
	basic_smatrix(T _xx, T _xy, T _yy) :
		xx(_xx), xy(_xy), yy(_yy)
		{ }
	template <class U> 
	basic_smatrix(const basic_smatrix<U> & v) :
		xx(T(v.xx)), xy(T(v.xy)), yy(T(v.yy))
		{ }
	/// convert from ordinary matrix by averaging its not-diagonal components
	template <class U> 
	explicit basic_smatrix(const basic_matrix<U> & v) :
		xx(T(v.xx)), xy(T(0.5) * T(v.xy + v.yx)), yy(T(v.yy))
		{ }
	basic_smatrix & identity(T diagval = 1) 
		{ xx = yy = diagval; xy = 0; return * this; }
	basic_smatrix & zero() 
		{ xx = yy = xy = 0; return * this; } 
public:
	basic_vector<T> operator ()(const basic_vector<T> & v) const {
		return basic_vector<T>(xx * v.x + xy * v.y,	xy * v.x + yy * v.y);
	}
public:
	T det() const
		{ return xx * yy - xy * xy; };
	operator bool() const ///< true if not degenerate
		{ return det() != 0; };
	basic_smatrix<T> inverse() const;
	/// uses precomputed value of determinant
	basic_smatrix<T> inverse(T det) const;
	basic_vector<T> inverse(const basic_vector<T> & v) const; ///< \return product of inverse matrix and given vector
	T tr() const ///< trace of the matrix
		{ return xx + yy; }
	T eigenvalue(int n) const; ///< n = 0 or 1
	basic_vector<T> eigenvector(int n, bool normalize = false, T * eigenvalue = 0) const; ///< n = 0 or 1
public:
    basic_smatrix & operator += (const basic_smatrix<T> & v)
		{ xx += v.xx; xy += v.xy; yy += v.yy; return * this; }
    basic_smatrix & operator += (T a) ///< adds identity matrix multiplied on a
		{ xx += a; yy += a; return * this; }
    basic_smatrix & operator -= (const basic_smatrix<T> & v)
		{ xx -= v.xx; xy -= v.xy; yy -= v.yy; return * this; }
    basic_smatrix & operator -= (T a)
		{ xx -= a; yy -= a; return * this; }
    basic_smatrix & operator *= (T a)
		{ xx *= a; xy *= a; yy *= a; return * this; }
    basic_smatrix & operator /= (T a)
		{ xx /= a; xy /= a; yy /= a; return * this; }
	bool operator == (const basic_smatrix<T> & v) const
		{ return xx == v.xx && xy == v.xy && yy == v.yy; }
	bool operator != (const basic_smatrix<T> & v) const
		{ return !(*this == v); }
};

template <class T>
inline basic_smatrix<T> basic_smatrix<T>::inverse(T det) const {
	return basic_smatrix<T>(yy, -xy, xx) / det;
}

template <class T>
inline basic_smatrix<T> basic_smatrix<T>::inverse() const {
	return inverse(det());
}

template <typename T>
inline basic_vector<T> basic_smatrix<T>::inverse(const basic_vector<T> & v) const {
	return basic_vector<T>(
		 yy * v.x - xy * v.y,
		-xy * v.x + xx * v.y
	) / det();
}

template <class T>
inline T basic_smatrix<T>::eigenvalue(int n) const {
	T a = 0.5 * (xx + yy);
	T b = sqrt(0.25 * sqr(xx - yy) + sqr(xy));
	return (n == 0) ? a - b : a + b;
}

template <class T>
inline basic_vector<T> basic_smatrix<T>::eigenvector(int n, bool normalize, T * eigenvalue) const {
	T a = 0.5 * (xx - yy);
	T b = sqrt(a*a + xy*xy);
	T c = a + b;
	T d = a - b;
	basic_vector<T> res;
	if (c*c > d*d)
		res = (n == 0) ? basic_vector<T>(xy, -c) : basic_vector<T>(c, xy);
	else
		res = (n == 0) ? basic_vector<T>(d, xy) : basic_vector<T>(xy, -d);
	if (normalize)
		res *= T(1.0 / absolute(res));
	if (eigenvalue)
		*eigenvalue = (n == 0) ? d + yy : c + yy;
	return res;
}

template <class T>
inline basic_smatrix<T> operator + (const basic_smatrix<T> & a, const basic_smatrix<T> & b) {
	basic_smatrix<T> t(a);
    return t += b;
}

template <class T>
inline basic_smatrix<T> operator + (const basic_smatrix<T> & v, T a) {
	basic_smatrix<T> t(v);
    return t += a;
}

template <class T>
inline basic_smatrix<T> operator - (const basic_smatrix<T> & a, const basic_smatrix<T> & b) {
	basic_smatrix<T> t(a);
    return t -= b;
}

//unary minus
template <class T>
inline basic_smatrix<T> operator - (const basic_smatrix<T> & a) {
	return basic_smatrix<T>(-a.xx, -a.xy, -a.yy);
}

template <class T>
inline basic_smatrix<T> operator - (const basic_smatrix<T> & v, T a) {
	basic_smatrix<T> t(v);
    return t -= a;
}

template <class T>
inline basic_smatrix<T> operator * (const basic_smatrix<T> & v, T a) {
	basic_smatrix<T> t(v);
    return t *= a;
}

template <class T>
inline basic_smatrix<T> operator * (T a, const basic_smatrix<T> & v) {
	basic_smatrix<T> t(v);
    return t *= a;
}

template <class T>
inline basic_vector<T> operator * (const basic_smatrix<T> & m, const basic_vector<T> & v) {
    return m(v);
}

//for commutative matrices only (product for general matrices see in g2/matrix.h")
template <typename T>
inline basic_smatrix<T> comm_product(const basic_smatrix<T> & a, const basic_smatrix<T> & b) {
    return basic_smatrix<T> (
		a.xx * b.xx + a.xy * b.xy,
		a.xx * b.xy + a.xy * b.yy,
		a.xy * b.xy + a.yy * b.yy
	);
}

template <class T>
inline basic_smatrix<T> operator / (const basic_smatrix<T> & v, T a) {
	basic_smatrix<T> t(v);
    return t /= a;
}

///tensor square (a*a)
template <class T>
inline basic_smatrix<T> tsquare(const basic_vector<T> & a) {
	return basic_smatrix<T>(
		a.x * a.x, //xx
		a.x * a.y, //xy
		a.y * a.y //yy
	);
}

///symmetric sum of tensor products (a*b + b*a)
template <class T>
inline basic_smatrix<T> sproduct(const basic_vector<T> & a, const basic_vector<T> & b) {
	return basic_smatrix<T>(
		2 * a.x * b.x, //xx
		a.x * b.y + a.y * b.x, //xy
		2 * a.y * b.y //yy
	);
}

///double convolution
template <class T>
inline T dscalar_product(const basic_smatrix<T> & v1, const basic_smatrix<T> & v2) 
	{ return v1.xx * v2.xx + v1.yy * v2.yy + 2 * v1.xy * v2.xy; }

// returns s * e * s 
template <class T>
inline basic_smatrix<T> duplex_product(const basic_smatrix<T> & e, const basic_smatrix<T> & s) 
{
	T sxxxx = s.xx * s.xx;
	T sxxxy = s.xx * s.xy;
	T sxxyy = s.xx * s.yy;
	T sxyxy = s.xy * s.xy;
	T sxyyy = s.xy * s.yy;
	T syyyy = s.yy * s.yy;
	return basic_smatrix<T>(
		sxxxx * e.xx + 2 * sxxxy * e.xy + sxyxy * e.yy, //xx
		sxxxy * e.xx + (sxyxy + sxxyy) * e.xy + sxyyy * e.yy, //xy
		sxyxy * e.xx + 2 * sxyyy * e.xy + syyyy * e.yy //yy
	);
}

} //namespace g2

template <class T>
class squarer<g2::basic_smatrix<T> > {
	T s;
public:
	explicit squarer(const g2::basic_smatrix<T> & v) : s(g2::dscalar_product(v, v)) { }
	operator T () const { return s; }
};

template <class T>
class absoluter<g2::basic_smatrix<T> > {
	T a;
public:
	explicit absoluter(const g2::basic_smatrix<T> & v) : a(::sqrt(sqr(v))) { }
	operator T () const { return a; }
};

#endif // glib_g2_smatrix_h
